library(readxl)
library(reshape2)
library(doBy)
library(dplyr)
library(stringr)
library(ggplot2)
library(ggforce)
library(vegan)
library(ggdendro)
library(patchwork)
library(ggrepel)
library(rjson)
library(car)
library(plyr)

ggplot theme setting

large <- theme(legend.title = element_text(size=15),
        legend.text = element_text(size=15),
        axis.title = element_text(size=15),
        axis.text = element_text(size=15),
        strip.text = element_text(size=15))

rotate <- theme(axis.text.x = element_text(angle = 90, hjust=1))

no_strip <- theme(strip.background = element_rect(colour=NA, fill=NA),
                  strip.text = element_text(colour=NA))

Preparing tanaid data

ta0 <- read_excel("../Data/Tanaidacea community.xlsx", sheet = 3)
ta <- ta0[, 1:2]
# Create factor "Microhabitat"
ta$Microhabitat <- str_replace(ta0$Microhabitat, " \\s*\\([^\\)]+\\)", "")
# Create factor "Replicate"
ta$Replicate <- str_extract(ta0$Microhabitat, "(?<=\\()\\d+(?=\\))")
ta <- cbind(ta, ta0[,-1:-3])
# Remove zero rows
ta <- ta[rowSums(ta0[,-1:-3])!=0, ]
ta$Site <- factor(ta$Site, levels=c("Shitiping", "Jihuei", "Jialulan"), labels = c("STP", "JH", "JLL"))
ta$Season <- factor(ta$Season, levels = c("SP", "SU"))

names(ta)[-1:-4] <- c("P. pangcahi", "Cyclopoapseudes sp.", "P. tagopilosus", "S. hansmuelleri", "I. multituberculata", "C. taitungensis", "P. setosa", "A. lenoprimorum","A. pedecerritulus", "T. nuwalianensis", "Z. shitipingensis", "Z. zorro")

ta$Microhabitat <- factor(ta$Microhabitat, 
                          levels = c("Amansia glomerata", "Asparagopsis sp.", "Asparagopsis taxiformis", "Caulerpa sp.", "Chlorodesmis sp.", "Corallina pilulifera", "Gelidium sp.", "Gravel, sand or silt", "Halimeda sp.", "Hypnea pannosa", "Hypnea sp.", "Jania sp.1", "Jania sp.2", "Mastophora rosea", "Padina sp.", "Plocamium sp.", "Tube of Eunice taoi"), 
                          labels = c("A. glomerata", "Asparagopsis sp.", "A. taxiformis", "Caulerpa sp.", "Chlorodesmis sp.", "C. pilulifera", "Gelidium sp.", "Gravel, sand, silt", "Halimeda sp.", "H. pannosa", "Hypnea sp.", "Jania sp.1", "Jania sp.2", "M. rosea", "Padina sp.", "Plocamium sp.", "Tube of E. taoi"))

# Only keep data for analysis
b <- ta[, -1:-4]
rownames(b) <- with(ta, paste(Site, Season, Microhabitat, Replicate, sep="_"))

Preparing microhabitat grain size data

cs <- read_excel("../Data/Grain-size composition.xlsx", sheet = 6) %>% as.data.frame

cs$Microhabitat <- factor(cs$Microhabitat, levels=c("Amansia glomerata", "Asparagopsis sp.", "Asparagopsis taxiformis", 
                                                    "Caulerpa sp.", "Chlorodesmis sp.", "Corallina pilulifera", "Gelidium sp.",
                                                    "sand", "Halimeda sp.", "Hypnea pannosa", "Jania sp.1", "Jania sp.2", 
                                                    "Mastophora rosea", "Padina sp.",  "polychaete tube"), 
                          labels = c("A. glomerata", "Asparagopsis sp.", "A. taxiformis", "Caulerpa sp.", "Chlorodesmis sp.", 
                                     "C. pilulifera", "Gelidium sp.", "Gravel, sand, silt", "Halimeda sp.", "H. pannosa",  
                                     "Jania sp.1", "Jania sp.2", "M. rosea", "Padina sp.", "Tube of E. taoi"))

# Average by site, season, and microhabitat
cs_avg <- aggregate(cs[,-1:-5], by=list(cs$Site, cs$Season, cs$Microhabitat), FUN=mean, na.rm=TRUE)
names(cs_avg)[1:3] <- c("Site", "Season", "Microhabitat")
cs_avg[is.na(cs_avg)] <- 0

cs[is.na(cs)] <- 0
cs_names <- with(cs,  paste(Site, Season, Microhabitat, Replicate, sep="_"))
ta_names <- with(ta,  paste(Site, Season, Microhabitat, Replicate, sep="_"))
env <- cs[match(ta_names, cs_names), -1:-5]
row.names(env) <- rownames(b)

# Logit transformation and normalized (divided by sd)
ens <- scale(logit(cs[, 6:9]/100))[match(ta_names, cs_names), ] %>% as.data.frame
row.names(ens) <- rownames(b)

Preparing CWB buoy data

cwb <-  fromJSON(file="../Data/C-B0050-001.json")

# Meta data 
meta <- lapply(cwb$cwbdata$Resources$Resource$Data$SeaSurfaceObs$Location, FUN=function(x)x[[1]])

# Monthly temperature
dat <- lapply(cwb$cwbdata$Resources$Resource$Data$SeaSurfaceObs$Location, FUN=function(x)x[[2]]$Monthly)

MonthlyTemperatureFun <- function(i){
  cbind(StationNameEN = lapply(meta, FUN=function(x)x$StationNameEN)[i] %>% as.character %>% as.factor,
        cbind(
          StationLatitude=lapply(meta, FUN=function(x)x$StationLatitude)[i] %>% as.numeric,
          StationLongitude=lapply(meta, FUN=function(x)x$StationLongitude)[i] %>% as.numeric,
          DataMonth=lapply(dat[[i]], FUN=function(x)x$DataMonth) %>% unlist %>% as.numeric,
          Maximum=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Maximum) %>% unlist %>% as.numeric,
          MaximumYear=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MaximumYear) %>% unlist %>% as.numeric,
          Mean=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Mean) %>% unlist %>% as.numeric,
          Minimum=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Minimum) %>% unlist %>% as.numeric,
          MinimumYear=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MinimumYear) %>% unlist %>% as.numeric,
          MinimumAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MinimumAnomaly) %>% unlist %>% as.numeric,
          MaximumAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MaximumAnomaly) %>% unlist %>% as.numeric,
          MeanAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MeanAnomaly) %>% unlist %>% as.numeric
          ) %>% as.data.frame
  )
}

# Hualien Buoy
hb <- MonthlyTemperatureFun(18) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")
# Chenggong
cg <- MonthlyTemperatureFun(8) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")
# Taitung Buoy
tb <- MonthlyTemperatureFun(34) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")

tmp <- rbind(hb, cg, tb)
tmp$Site <- factor(tmp$StationNameEN, levels=c("Hualien Buoy", "Chenggong", "Taitung Buoy"), labels = c("STP", "JH", "JLL"))
tmp$Season <- factor(tmp$DataMonth, levels=c(4, 8), labels=c("SP", "SU"))

tmp_names <- with(tmp, paste(Site, Season, sep="_"))
ta_names <- with(ta,  paste(Site, Season, sep="_"))

env$Temperature <- tmp$Mean[match(ta_names, tmp_names)]
# Log 10 transformation and normalize (divided by sd)
ens$Temperature <- scale(log10(env$Temperature))

By replicates

Cluster Analysis

# Hierarchical Clustering based on square root transformed abundance data 
# Agglomeration used group average (= UPGMA) method
d <- vegdist(b^0.25)
hc <- hclust(d, method="average")
#plot(x = hc, labels =  row.names(hc), cex = 0.5, hang=-1)
#rect.hclust(tree = hc, k = 4)
fac <- ta[, 1:4]
fac$Group <- cutree(hc, 4)
# Reorder group from north to south as Groups 1, 2, 3
fac$Group <- factor(fac$Group, levels=c(1,2,3,4), labels = c(2, 3, 1, 4))

# Convert dendrogram to ggplot style
dhc <- as.dendrogram(hc)
ghc    <- dendro_data(dhc, type="rectangle") 

# Merge sample info to dendrogram label data frame
ord <- match(label(ghc)$label, rownames(b))
ghc[["labels"]]   <- cbind(ghc[["labels"]], fac[ord,])
# Group membership

# Group 1 (Shitiping = 82.6%)
table(subset(fac, Group==1)$Site)
## 
## STP  JH JLL 
##  19   2   2
# Group 2 (Jihuei = 52.6%, Jialulan = 28.9%)
table(subset(fac, Group==2)$Site)
## 
## STP  JH JLL 
##   7  20  11
# Group 3 (Jialulan = 77.6%, Jihuei = 17.2%)
table(subset(fac, Group==3)$Site)
## 
## STP  JH JLL 
##   3  10  45

Nonmetric Multidimensional Scaling

# Remove outlier JH_AP_Jania sp.2_1
md <- metaMDS(vegdist(b[-23,]^0.25))
## Run 0 stress 0.1395381 
## Run 1 stress 0.1493668 
## Run 2 stress 0.1431431 
## Run 3 stress 0.1454937 
## Run 4 stress 0.142695 
## Run 5 stress 0.1418538 
## Run 6 stress 0.1418331 
## Run 7 stress 0.1380858 
## ... New best solution
## ... Procrustes: rmse 0.04631679  max resid 0.1232183 
## Run 8 stress 0.1383089 
## ... Procrustes: rmse 0.01505367  max resid 0.07531758 
## Run 9 stress 0.1397265 
## Run 10 stress 0.14193 
## Run 11 stress 0.1406395 
## Run 12 stress 0.1397602 
## Run 13 stress 0.1382344 
## ... Procrustes: rmse 0.004775539  max resid 0.04929032 
## Run 14 stress 0.1383071 
## ... Procrustes: rmse 0.01515002  max resid 0.07518309 
## Run 15 stress 0.1398386 
## Run 16 stress 0.1380858 
## ... New best solution
## ... Procrustes: rmse 3.250285e-05  max resid 0.0002187851 
## ... Similar to previous best
## Run 17 stress 0.1395381 
## Run 18 stress 0.1414666 
## Run 19 stress 0.1389654 
## Run 20 stress 0.1425775 
## *** Best solution repeated 1 times
stress <- paste("Stress = ", deparse(round(md$stress,2)))

keep <- names(colSums(b[-23,])%>%sort(decreasing = TRUE))[1:5]

blab <- wascores(md$points, b[-23,]^0.25, expand = TRUE)[keep,] %>% as.data.frame
blab$label <- paste("italic(", sub(" ", "~~", rownames(blab)), ")", sep="")

# Only fits environmental vector to the replicate with grain size data
keep <- !env[-23, 1:4]%>%rowSums%>%is.na
set.seed(100)
(fit <- envfit(md$points[keep,], env[-23,][keep,]))
## 
## ***VECTORS
## 
##                       MDS1     MDS2     r2 Pr(>r)    
## Silt2Clay          0.29515  0.95545 0.0264  0.475    
## Fine2VeryFine      0.53410  0.84542 0.1250  0.031 *  
## Medium             1.00000  0.00267 0.2581  0.001 ***
## Coarse2VeryCoarse  0.88051 -0.47403 0.0912  0.078 .  
## Temperature       -0.69724  0.71683 0.0690  0.158    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
evec <- as.data.frame(scores(fit, display = "vectors")) * ordiArrowMul(fit, fill = 1.5)
evec$label <- row.names(evec)

# Mantel tests on each variables
man <- list()
man[[1]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Silt2Clay"]))
man[[2]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Fine2VeryFine"]))
man[[3]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Medium"]))
man[[4]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Coarse2VeryCoarse"]))
man[[5]] <- mantel(vegdist(b[-23,]^0.25), dist(env[-23, "Temperature"]))

evec$mantel_r <- lapply(man, FUN=function(x)x$statistic) %>% unlist
evec$signif <- lapply(man, FUN=function(x)x$signif) %>% unlist
#evec$pvals <- cut(evec$signif, breaks=c(0, 0.01, 0.05, 0.1, 1), labels=c("<0.01", "<0.05", "<0.1", "<1")) %>% factor
evec$pvals <- factor(evec$signif)

print(evec)
##                         MDS1         MDS2             label    mantel_r signif
## Silt2Clay          0.1415600  0.458253191         Silt2Clay -0.05470974  0.827
## Fine2VeryFine      0.5575815  0.882594237     Fine2VeryFine  0.07029091  0.083
## Medium             1.5000000  0.003998335            Medium  0.21036458  0.001
## Coarse2VeryCoarse  0.7850268 -0.422627016 Coarse2VeryCoarse  0.11119698  0.006
## Temperature       -0.5406990  0.555889465       Temperature  0.13973555  0.001
##                   pvals
## Silt2Clay         0.827
## Fine2VeryFine     0.083
## Medium            0.001
## Coarse2VeryCoarse 0.006
## Temperature       0.001
# Merge MDS to environmental data frame
p2 <- ggplot(data=cbind(md$points, fac[-23,]), aes(x=MDS1, y=MDS2))+
  geom_point(alpha=1, size=5, stroke=1, aes(colour=Site, shape=Season))+
  geom_mark_hull(concavity = 10, aes(group=Group), colour="black", linetype=2)+
  geom_segment(data= evec, aes(x=0, y=0, xend=MDS1, yend=MDS2, linewidth=pvals, linetype=pvals),
               arrow = arrow(length=unit(.4, 'cm')), colour="red")+
  annotate("text", x=-0.9, y=1.1, label=stress, size=5) +
  annotate("text", x=1, y=0, label="1", size=10) +
  annotate("text", x=-0.5, y=-0.6, label="2", size=10) +
  annotate("text", x=-0.2, y=0.5, label="3", size=10) +
  annotate("text", x=1.7, y=1.2, label="b", size=7) +
  scale_shape_manual(values=c(1, 19))+
  scale_color_viridis_d()+
  scale_linewidth_manual(values=c(2, 1, 0.5, 0.5))+
  scale_linetype_manual(values=c(1, 1, 1, 2))+
  geom_label_repel(data=blab, aes(x=MDS1, y=MDS2, label=label), colour="black", 
             fill=gray(1, 0.6), size=4, label.padding = unit(0.25, "lines"), parse=T)+
  geom_label_repel(data=evec, aes(x=MDS1*1.2, y=MDS2*1.2, label=label), colour="red", 
             fill=gray(1, 0.6), size=4, label.padding = unit(0.25, "lines"), force=2)+
  labs(linewidth="p-value", linetype="p-value")+
  theme_bw() %+replace% large
ggdendrogram(hc, rotate = FALSE)+
  geom_point(data=label(ghc), aes(x, y, colour=Site, shape=Season), size=2, stroke=0.5)+
  geom_hline(yintercept = 0.8, linetype=2, colour="red")+
  annotate("text", x = 95, y = 0.85, label = "2", size = 5)+
  annotate("text", x = 42, y = 0.85, label = "3", size = 5)+
  annotate("text", x = 10, y = 0.85, label = "1", size = 5)+
  scale_shape_manual(values=c(1, 19))+
  scale_color_viridis_d()+
  labs(x= "Replicate",y = "Bray-Curtis Dissimilarity")+
  theme_bw() %+replace% theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.line.y = element_line(colour="black")) %+replace%  rotate

m <- acast(dat=ta, Microhabitat~Site+Season, fun.aggregate = length) 
hc <- vegdist (m) %>% hclust(method="average")
clust_names <- row.names(m)[rev(hc$order)]
out <- cbind(md$points, fac[-23,])
out$Microhabitat <- factor(out$Microhabitat, levels=clust_names)

# Color code the microhabitat
cols <- c("#222222", "#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", rep("#c2b280", 7), rep("#848482", 3),"#008856")
shapes <- c(rep(19, 6), 0:6, 15, 17, 18:19)


# MDS for each site and season
md <- lapply(splitBy(~Site+Season, ta), FUN=function(x)metaMDS(x[,-1:-4]))
## Wisconsin double standardization
## Run 0 stress 0.09177207 
## Run 1 stress 0.09177187 
## ... New best solution
## ... Procrustes: rmse 0.000338808  max resid 0.001136196 
## ... Similar to previous best
## Run 2 stress 0.09240116 
## Run 3 stress 0.1102589 
## Run 4 stress 0.09185124 
## ... Procrustes: rmse 0.07012022  max resid 0.1591274 
## Run 5 stress 0.09901149 
## Run 6 stress 0.09185037 
## ... Procrustes: rmse 0.06994872  max resid 0.1586776 
## Run 7 stress 0.1204544 
## Run 8 stress 0.09857996 
## Run 9 stress 0.09058073 
## ... New best solution
## ... Procrustes: rmse 0.1173281  max resid 0.306095 
## Run 10 stress 0.09240098 
## Run 11 stress 0.1198589 
## Run 12 stress 0.09185052 
## Run 13 stress 0.0955431 
## Run 14 stress 0.09050981 
## ... New best solution
## ... Procrustes: rmse 0.08364446  max resid 0.2192915 
## Run 15 stress 0.0909706 
## ... Procrustes: rmse 0.0747375  max resid 0.2755401 
## Run 16 stress 0.1253892 
## Run 17 stress 0.125494 
## Run 18 stress 0.09050979 
## ... New best solution
## ... Procrustes: rmse 0.0006690435  max resid 0.001675329 
## ... Similar to previous best
## Run 19 stress 0.09058075 
## ... Procrustes: rmse 0.08327585  max resid 0.2200089 
## Run 20 stress 0.09058272 
## ... Procrustes: rmse 0.08290963  max resid 0.2208044 
## *** Best solution repeated 1 times
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1345332 
## Run 1 stress 0.1480562 
## Run 2 stress 0.1468133 
## Run 3 stress 0.1551481 
## Run 4 stress 0.1383558 
## Run 5 stress 0.153372 
## Run 6 stress 0.1322342 
## ... New best solution
## ... Procrustes: rmse 0.09070836  max resid 0.298377 
## Run 7 stress 0.1403992 
## Run 8 stress 0.1755181 
## Run 9 stress 0.1678867 
## Run 10 stress 0.1564364 
## Run 11 stress 0.1442084 
## Run 12 stress 0.1524216 
## Run 13 stress 0.1361397 
## Run 14 stress 0.1403993 
## Run 15 stress 0.17642 
## Run 16 stress 0.153372 
## Run 17 stress 0.1346216 
## Run 18 stress 0.1678868 
## Run 19 stress 0.1550996 
## Run 20 stress 0.139706 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1208441 
## Run 1 stress 0.09427342 
## ... New best solution
## ... Procrustes: rmse 0.1234327  max resid 0.2720241 
## Run 2 stress 0.1098328 
## Run 3 stress 0.09427323 
## ... New best solution
## ... Procrustes: rmse 7.529209e-05  max resid 0.0002953634 
## ... Similar to previous best
## Run 4 stress 0.1382285 
## Run 5 stress 0.1349709 
## Run 6 stress 0.09332658 
## ... New best solution
## ... Procrustes: rmse 0.03245163  max resid 0.1489107 
## Run 7 stress 0.1196955 
## Run 8 stress 0.09427334 
## Run 9 stress 0.1212999 
## Run 10 stress 0.125197 
## Run 11 stress 0.1175232 
## Run 12 stress 0.1265351 
## Run 13 stress 0.1170811 
## Run 14 stress 0.09332658 
## ... New best solution
## ... Procrustes: rmse 0.000273151  max resid 0.0009519536 
## ... Similar to previous best
## Run 15 stress 0.1099111 
## Run 16 stress 0.09427323 
## Run 17 stress 0.09332666 
## ... Procrustes: rmse 0.0003160555  max resid 0.001152594 
## ... Similar to previous best
## Run 18 stress 0.1145061 
## Run 19 stress 0.1054994 
## Run 20 stress 0.1294716 
## *** Best solution repeated 2 times
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.138881 
## Run 1 stress 0.1380855 
## ... New best solution
## ... Procrustes: rmse 0.02290013  max resid 0.09354521 
## Run 2 stress 0.2044978 
## Run 3 stress 0.2363496 
## Run 4 stress 0.1750268 
## Run 5 stress 0.2121027 
## Run 6 stress 0.2069126 
## Run 7 stress 0.1518363 
## Run 8 stress 0.2125752 
## Run 9 stress 0.1393694 
## Run 10 stress 0.216466 
## Run 11 stress 0.1553091 
## Run 12 stress 0.133009 
## ... New best solution
## ... Procrustes: rmse 0.08724517  max resid 0.3626722 
## Run 13 stress 0.1638621 
## Run 14 stress 0.1400827 
## Run 15 stress 0.2068781 
## Run 16 stress 0.2026245 
## Run 17 stress 0.1564282 
## Run 18 stress 0.1526192 
## Run 19 stress 0.2163723 
## Run 20 stress 0.2174266 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     15: stress ratio > sratmax
##      5: scale factor of the gradient < sfgrmin
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.05720906 
## Run 1 stress 0.0572091 
## ... Procrustes: rmse 0.0002015069  max resid 0.0003602265 
## ... Similar to previous best
## Run 2 stress 0.0572092 
## ... Procrustes: rmse 0.0002634092  max resid 0.0004780251 
## ... Similar to previous best
## Run 3 stress 0.05720917 
## ... Procrustes: rmse 0.0002193702  max resid 0.0004674462 
## ... Similar to previous best
## Run 4 stress 0.05720909 
## ... Procrustes: rmse 7.852231e-05  max resid 0.00014762 
## ... Similar to previous best
## Run 5 stress 0.05720917 
## ... Procrustes: rmse 0.0002701912  max resid 0.0004919481 
## ... Similar to previous best
## Run 6 stress 0.05720914 
## ... Procrustes: rmse 0.0002475757  max resid 0.0004502463 
## ... Similar to previous best
## Run 7 stress 0.05720912 
## ... Procrustes: rmse 0.0002307228  max resid 0.0004188804 
## ... Similar to previous best
## Run 8 stress 0.08766255 
## Run 9 stress 0.08093205 
## Run 10 stress 0.08093206 
## Run 11 stress 0.1086705 
## Run 12 stress 0.05720904 
## ... New best solution
## ... Procrustes: rmse 4.555359e-05  max resid 7.944265e-05 
## ... Similar to previous best
## Run 13 stress 0.05720907 
## ... Procrustes: rmse 0.0001330564  max resid 0.0002378443 
## ... Similar to previous best
## Run 14 stress 0.08093209 
## Run 15 stress 0.05720903 
## ... New best solution
## ... Procrustes: rmse 3.619008e-05  max resid 0.0001015064 
## ... Similar to previous best
## Run 16 stress 0.05720904 
## ... Procrustes: rmse 2.007738e-05  max resid 5.202092e-05 
## ... Similar to previous best
## Run 17 stress 0.07810263 
## Run 18 stress 0.08093205 
## Run 19 stress 0.08093206 
## Run 20 stress 0.05720908 
## ... Procrustes: rmse 0.0001108149  max resid 0.0002103537 
## ... Similar to previous best
## *** Best solution repeated 3 times
## Wisconsin double standardization
## Run 0 stress 6.522281e-33 
## Run 1 stress 0 
## ... New best solution
## ... Procrustes: rmse 0.1307195  max resid 0.2126991 
## Run 2 stress 0 
## ... Procrustes: rmse 0.1368719  max resid 0.2230256 
## Run 3 stress 0 
## ... Procrustes: rmse 0.1450721  max resid 0.287975 
## Run 4 stress 0 
## ... Procrustes: rmse 0.1626683  max resid 0.2638379 
## Run 5 stress 0 
## ... Procrustes: rmse 0.09148311  max resid 0.1643084 
## Run 6 stress 1.910875e-05 
## ... Procrustes: rmse 0.149761  max resid 0.2413231 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1231937  max resid 0.2250026 
## Run 8 stress 0 
## ... Procrustes: rmse 0.1865359  max resid 0.3227937 
## Run 9 stress 0 
## ... Procrustes: rmse 0.1331486  max resid 0.2704025 
## Run 10 stress 0 
## ... Procrustes: rmse 0.115681  max resid 0.1553413 
## Run 11 stress 0 
## ... Procrustes: rmse 0.1302713  max resid 0.2177031 
## Run 12 stress 6.264849e-05 
## ... Procrustes: rmse 0.1032593  max resid 0.1609292 
## Run 13 stress 0 
## ... Procrustes: rmse 0.139783  max resid 0.2708793 
## Run 14 stress 0 
## ... Procrustes: rmse 0.1839108  max resid 0.3134858 
## Run 15 stress 0 
## ... Procrustes: rmse 0.1270628  max resid 0.1964893 
## Run 16 stress 0 
## ... Procrustes: rmse 0.172901  max resid 0.2816818 
## Run 17 stress 0 
## ... Procrustes: rmse 0.1588443  max resid 0.2648712 
## Run 18 stress 0 
## ... Procrustes: rmse 0.1501897  max resid 0.2424697 
## Run 19 stress 0 
## ... Procrustes: rmse 0.1613569  max resid 0.2847132 
## Run 20 stress 0 
## ... Procrustes: rmse 0.1396997  max resid 0.2310131 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin
# Stress levels
st <- lapply(md, FUN=function(x)x$stress) %>% ldply
st <- cbind(strsplit(st$.id, split="[|]") %>% ldply, st[, -1])
names(st) <- c("Site", "Season", "stress")
st$stress <- paste("Stress = ", round(st$stress,2))
st$Site <- factor(st$Site, levels=c("STP", "JH", "JLL"))
# MDS scores
md <- cbind(ta[,1:4], lapply(md, FUN=function(x)x$points) %>% ldply)
md$Microhabitat <- factor(md$Microhabitat, levels = clust_names)

ggplot(data=md, aes(x=MDS1, y=MDS2, colour=Microhabitat, shape=Microhabitat))+
  geom_point(size=3)+
  geom_text(data = st, inherit.aes = FALSE, aes(x=Inf, y = Inf, label = stress), vjust=1.2, hjust=1.05)+
  scale_colour_manual(values=cols)+
  scale_shape_manual(values=shapes)+
  facet_grid(Site~Season, scales="free")+
  theme_bw() %+replace% large

Aggregate by site-season

ta3 <- aggregate(b, by = list(ta$Site, ta$Season), FUN=mean)
names(ta3)[1:2] <- c("Site", "Season")
b3 <- ta3[, -1:-2]
rownames(b3) <- with(ta3, paste(Site, Season, sep="_"))

# Average by site and season
cs_avg2 <- aggregate(cs[,-1:-5], by=list(cs$Site, cs$Season), FUN=mean, na.rm=TRUE)
names(cs_avg2)[1:2] <- c("Site", "Season")
cs_avg2
##   Site Season  Silt2Clay Fine2VeryFine    Medium Coarse2VeryCoarse
## 1   JH     SP 0.54811731      8.016581 10.550036          55.88527
## 2  JLL     SP 0.72627053     12.150971  3.015219          12.67897
## 3  STP     SP 0.07468493      9.590777 12.213665          42.40659
## 4   JH     SU 0.93316520     10.156021  9.265062          59.64575
## 5  JLL     SU 0.65368487     15.726078  5.927816          33.94242
## 6  STP     SU 2.85256410     16.025641 18.878205          62.24359

Cluster Analysis

# Hierarchical Clustering based on square root transformed abundance data 
# Agglomeration used group average (= UPGMA) method
(d <- vegdist(b3^0.25))
##           STP_SP     JH_SP    JLL_SP    STP_SU     JH_SU
## JH_SP  0.3183515                                        
## JLL_SP 0.3953006 0.3404144                              
## STP_SU 0.3071157 0.3873459 0.4548924                    
## JH_SU  0.3853953 0.1707948 0.2189590 0.4929781          
## JLL_SU 0.3759794 0.1824998 0.2880368 0.3637466 0.1863058
hc <- hclust(d, method="average")
fac <- ta3[, 1:2]

# Convert dendrogram to ggplot style
dhc <- as.dendrogram(hc)
ghc    <- dendro_data(dhc, type="rectangle") 

# Merge sample info to dendrogram label data frame
ord <- match(label(ghc)$label, rownames(b3))
ghc[["labels"]]   <- cbind(ghc[["labels"]], fac[ord,])

p5 <- ggdendrogram(hc, rotate = FALSE)+
  geom_point(data=label(ghc), aes(x, y, colour=Site, shape=Season), size=4, stroke=1.3)+
  annotate("text", x=6, y=0.4, label="a", size=7) +
  scale_shape_manual(values=c(1, 19))+
  scale_color_viridis_d()+
  labs(x= "",y = "Bray-Curtis Dissimilarity")+
  theme_bw() %+replace% theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(),axis.line.y = element_line(colour="black")) %+replace% large

Comibining site-season clustering with replicate-level nMDS

p5/p2+plot_layout(nrow = 2, heights = c(1, 3))